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Abstract 

Theory considers the covariation of seasonal life-history traits as an optimal re- 
action norm, implying that deviating from this reaction norm reduces fitness. 
However, the estimation of reaction-norm properties (i.e., elevation, linear slope, 
and higher order slope terms) and the selection on these is statistically challeng- 
ing. We here advocate the use of random regression mixed models to estimate 
reaction-norm properties and the use of bivariate random regression to estimate 
selection on these properties within a single model. We illustrate the approach by 
random regression mixed models on 1115 observations of clutch sizes and lay- 
ing dates of 361 female Ural owl Strix uralensis collected over 31 years to show 
that (1) there is variation across individuals in the slope of their clutch size-laying 
date relationship, and that (2) there is selection on the slope of the reaction norm 
between these two traits. Hence, natural selection potentially drives the negative 
covariance in clutch size and laying date in this species. The random-regression 
approach is hampered by inability to estimate nonlinear selection, but avoids a 
number of disadvantages (stats-on-stats, connecting reaction-norm properties to 
fitness) . The approach is of value in describing and studying selection on behavioral 
reaction norms (behavioral syndromes) or life-history reaction norms. The ap- 
proach can also be extended to consider the genetic underpinning of reaction-norm 
properties. 



Introduction 

Reaction-norm theory states that covariance of life-history 
traits occurs when individuals experience a heterogeneous 
environment, where the optimal combination of life-history 
traits depends on the environmental context an individual 
finds itself in (Stearns 1992; Kozlowski 1992). For example, 
in practically all temperate avian species, birds that lay later 
in the breeding season have a lower clutch size than individ- 
uals which produce earlier in the breeding season (Klomp 
1970). The negative covariance of these two life-history traits 
(clutch size and laying date) is generally viewed as a reac- 
tion norm in response to environmental heterogeneity. Under 
poor environmental conditions, it is optimal for an individ- 
ual's fitness to reproduce late in the season and produce a 
small number of eggs compared to individuals that experi- 
ence good conditions (Fig. 1; Daan et al. 1990; Rowe et al. 
1994; Brommer et al. 2002a). Consider, for example, a lin- 
ear increase over time t in the potential clutch size C of 



individual i as 

C,(f) = C„, ; +at, (1) 

where potential clutch size C can be considered a common 
currency to express differences across individuals in the en- 
vironmental conditions they experience, for example, due to 
the food supply in their territory. All interindividual varia- 
tion in this example is fully due to variation in the initial 
potential clutch size experienced by individual i, Co,, (but see 
Rowe et al. 1994 and Brommer et al. 2002a for other scenar- 
ios). The reproductive value V of an egg produced at time t 
declines as 

V(t)=V 0 -/3t. (2) 

The optimal clutch size C* and laying date t * are then found 
by maximizing the product C(f) V(t) for every individual 
which generalizes (details in Brommer et al. 2002a) to 

C*(f) = aV 0 //3 - at*. (3) 
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Selection on Reaction-Norm Properties 

Hence, the individually optimized reaction norm is a de- 
cline in clutch size with advancing laying date at rate a. Al- 
though this concept is here illustrated with clutch size and 
laying date, the same arguments apply to other seasonal life- 
history traits such as, for example, size at maturation (re- 
lated to fecundity) and seasonal developmental time in in- 
sects (Kawecki and Stearns 1993). 

The above example of the theory behind reaction norms 
implies that individuals that do not follow the optimal re- 
action norm suffer reduced fitness compared to those that 
follow the norm. In the above example, if an individual ad- 
vances laying date relative to the reaction norm (gray squares 
in Fig. 1), fecundity is reduced and this reduction exceeds 
the gain in reproductive value of the offspring produced. If 
the individual delays and produces a larger clutch size rela- 
tive to the reaction norm (gray dots in Fig. 1), fitness costs 
occur in terms of reduced reproductive value of the offspring 
produced which outweigh the increase in fecundity. Thus, 
evaluation of empirical reaction-norm data should reveal se- 
lection on the properties describing the reaction norm. In 
particular, theory predicts that the slope of the trait-trait re- 
lationship should be under selection. Despite the intuitive 
appeal of the theory, empirical demonstration that a certain 




Time 



Figure 1. Illustration of the main theoretical background of the reaction- 
norm concept in terms of two seasonal life-history traits. As the season 
advances (Time), environmental conditions and thus the potential clutch 
size C increase (Equation [1]; dotted line), but individuals experience dif- 
ferent trajectories characterized by the initial condition C 0 . Because the 
reproductive value of offspring V declines seasonally (Equation [2]), there 
is an optimal switching curve (solid line) that describes when an individ- 
ual should stop increasing in condition and reproduce (filled dots). This 
switching curve is the optimal reaction norm maximizing V{t)C(t), show- 
ing a seasonal decrease in reproduction (Equation [3]). Timing repro- 
duction earlier has fitness costs, because reproductive output decreases 
(gray squares), more so than offspring reproduction increases. Delaying 
reproduction increases reproductive output (gray dots), but has fitness 
costs because late-produced offspring are of lower value. 
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covariation of traits is under selection is challenging. This is 
because it requires two steps: 

( 1 ) Quantification of the variation across individuals or geno- 
types in their reaction-norm properties. In its simplest 
form, a reaction norm is a linear function (Trait Y = 
elevation + slope x Trait X), thereby requiring estima- 
tion of the variance of at least two properties. Whereas 
the "elevation" is clearly related to the mean value of the 
trait, the "slope" (within-individual covariance between 
two traits) is an emergent property of the observed data 
and is therefore sensitive to choice of statistical procedure 
used to estimate it. Some of these procedures will have 
more power to describe variation in slope than others. 

(2) Selection on a reaction norm requires the calculation of 
the fitness of these reaction norm properties, rather than 
the fitness of the expression of a trait itself (for more dis- 
cussion on this notion, see Weis and Gorman (1990), and 
Spitze and Sadler (1996)). Because fitness is a trait of the 
individual, this implies that, in statistical terms, a descrip- 
tion of reaction-norm fitness must allow for the covari- 
ance between fitness and reaction norm slope, where the 
latter is an emergent property (a "meta-trait") defined by 
the within-individual covariance of two traits. 

Brommer et al. (2003) noted that for labile traits for which 
expression changes across an individual's lifetime, measure- 
ment of life-history traits (under varying environmental con- 
ditions) and lifetime fitness can be obtained for each individ- 
ual. Hence, selection on plasticity can be calculated for such 
traits. Nevertheless, calculation of selection on plasticity has 
thus far been based on regression of estimates of an individ- 
ual's reaction-norm slope (describing the covariance between 
two traits or between a trait and an environmental value) on 
an individual's fitness, thereby doing stats-on-stats, either 
on slope estimates derived from linear regression (Brommer 
et al. 2003; Nussey et al. 2005a) or derived from a linear mixed 
model (as Best Linear Unbiased Predictor [BLUP] values, 
Brommer et al. 2005; Nussey et al. 2005b). Linear regression 
of data on each individual requires a certain minimal num- 
ber of repeated records per individual (e.g., five in Brommer 
et al. 2005; four in Nussey et al. 2005a), thereby discard- 
ing a potentially sizeable proportion of the data consisting 
of observations on individuals with relatively few repeated 
records. In terms of selection analysis, this introduces a pos- 
sible bias as only the longest lived individuals are included. 
Mixed model BLUP values of a response variable are based 
on the assumption that all covariates that may have an effect 
on the response variable have been included (Pinheiro and 
Bates 2000). Hence, using BLUP values to explore covaria- 
tion with another trait (e.g., fitness), violates the assumption 
under which they were generated in the first place. Typically, 
both approaches ignore the often large uncertainty inherent 
in estimating slope when quantifying selection on the slope 
estimate, which may produce spurious estimates of selection. 
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Here, we use random regression mixed models to calculate 
both the reaction-norm properties and selection on these in 
an integrative manner that circumvents the above-described 
problems. Random regression is an implementation of the 
infinite-dimensional model (e.g., Kirkpatrick and Heckman 
1989), which presents a simplification of the character state 
approach by allowing trait variation to be described as a func- 
tion of a covariate. As all mixed models, random regression 
is robust for analyzing unbalanced data, which data from 
natural populations often are (Crawley 2002). Here, we ar- 
gue that a mixed model random regression can be used to 
explore variation in the properties underlying a trait-trait re- 
action norm. Random regression mixed models have mostly 
been used in the study of trait-environment relationships 
in natural populations, both on the phenotypic and genetic 
levels (reviewed by Nussey et al. 2007), but have — to our 
knowledge — not been implemented in the study of covari- 
ance between traits. Furthermore, we demonstrate how mul- 
tivariate random regression can be used to estimate selection 
on reaction-norm properties. Multivariate random regres- 
sion models have been rarely used in ecology, with the ex- 
ception of exploration on how the covariance between traits 
changes over an environmental gradient (e.g., Robinson et al. 
2008) and in formal comparison of trait-environment reac- 
tion norms in two populations (Husby et al. 2010). 

We illustrate our approach by a consideration of selection 
on the properties describing the clutch size-laying date re- 
lationship. Our study species, the Ural owl (Strix uralensis 
Pall.), experiences a highly variable environment, and shows 
high plasticity in clutch size (one to eight eggs) and seasonal 
timing of laying (more than 2 months) (e.g., Pietiainen 1989). 
Previous work (based on linear regression and analysis of co- 
variance) has shown that there is variation across Ural owl 
females in their clutch size-laying date relationship in terms 
of elevation and slope (Brommer et al. 2003). Here, we quan- 
tify the variation in these properties and selection on them 
using random-regression analysis. We discuss the applicabil- 
ity of this method in the study of phenotypic integration of 
behavioral and life-history traits. 



Methods: Random-Regression 
Approach 

We outline the model approach in terms of the analysis of 
clutch size and laying date. We here adhere to this specific 
scenario, instead of a more generic one, in order to facilitate 
understanding the approach and linking it to the example sce- 
nario presented in this paper and the code for implementing 
this model, which is presented in the Supporting Informa- 
tion. Nevertheless, the approach is applicable to also other 
combinations of traits that are expressed repeatedly during 
an individual's lifetime. 



Model construction and evaluation was performed in two 
stages. We started by modeling the clutch size-laying date 
relationships, such that 

C yr ,i = V-c+ AgeF ( + bd yrJ + yr + f(ind c x ■ d yrJ ) + E yrJ , 

(4) 

where c yr j denotes the clutch size of female i in year yr. The 
"/x c " fits the overall fixed effects mean clutch size, AgeF yr j 
is a fixed effect that denotes the factor age of individual i in 
year yr, d yr j the laying date of female i in year yr, and b the 
fixed-effect slope of clutch size as a function of laying date. 
Any annual variation that is not explained by the fixed effects 
is modeled by the random effect yr. The random regression 
element is given by the polynomial function f(ind c x ; , d yr j) 
in Equation (4), which specifies the individual-specific devi- 
ations from the fixed-effect slope b of clutch size on laying 
date. For each individual i, a polynomial function of increas- 
ing order x is specified. Thus, when x = 0, the variance 
across individuals (ind 0 ) is estimated at the average laying 
date value. When x = 1, the variance in the coefficients ind 0 
(elevation) and the variance in indi (slope) of the function 
indo + ind\d year is estimated and the covariance between 
these, and so on for higher order polynomials. Equation (4) 
is a standard random regression model, with the exception 
that typically an environmental variable is used as explana- 
tory variable (e.g., Schaeffer 2004; Nussey et al. 2007). 

To decide the order x of the random regression, we assumed 
that the most parsimonious random regression model was 
reached when higher order polynomials did not achieve a sig- 
nificant increase in log-likelihood. That is, the order x of the 
random regression polynomial function was increased step- 
wise, and its significance was tested by a likelihood ratio test, 
which is two times the difference in log-likelihood between 
hierarchically nested models, tested against a chi-square dis- 
tribution assuming that the degrees of freedom are given by 
the additional number of (co)variances estimated. Based on 
the most parsimonious random regression model, estimates 
of variance in clutch size at each laying date and its approxi- 
mate standard error can be calculated following Fischer et al. 
(2004). 

Selection on reaction-norm properties 

Having determined the most parsimonious order of the poly- 
nomial random regression, we extended the univariate ran- 
dom regression (Equation [4]) to a bivariate one by in addi- 
tion to Equation (4) considering that 

Wi = n w + CohortF, + f (ind^, d yrJ ) , (5) 

where w; is the relative fitness of individual The fi n , denotes 
the overall fixed effect mean fitness, and CohortF ; a fixed ef- 
fect denoting the year of first breeding of individual i, which 
is entered to correct for temporal variation in fitness and 
the possible truncation of lifetime fitness in the last cohorts 



©2012 The Authors. Published by Blackwell Publishing Ltd. 



697 



Selection on Reaction-Norm Properties 



J. E. Brommer ef al. 



considered in the analysis. The bivariate element in the ran- 
dom regression is specified by the polynomial f{ind c x d yr j). 
Thus, for order x, the variances in i nd™ and i nd c x are both es- 
timated, in addition to all possible covariances between these 
parameters. Nevertheless, because each individual has only 
one estimate of lifetime fitness, only the variance in elevation 
(indg) is estimable and all other variances (including the 
residual variance) must be constrained to zero. Thus, for a 
first-order random regression, the mixed model procedure 
needs to estimate the (co)variances (written in matrix form, 
omitting the upper diagonal that is symmetric to the lower 
diagonal). 

var(ind^) 
cov{ind c 0 , ind[) var{ind[) 
cov{ind c 0 , ind™) cov{ind[, ind™) var(indg) 

0 0 0 0_ 

(6) 

where cov(indg, ind™) denotes the selection differential 
on the elevation of the clutch size-laying date relationship 
(SeWation) and cov{ind[, ind™) the selection differential on 
the slope of the clutch size-laying date relationship (S s i ope ). 
These selection differentials can be expressed as a selection 
gradient by dividing them with variind^) and var{ind[), 
respectively, or as selection intensity (standardized selection 
gradient) by dividing them with the square root of these 
variances. 

Higher order polynomials can be introduced as an ex- 
tension of this matrix. Thus, bivariate random regression 
allows estimating directional selection on the properties of 
a polynomial individual-specific relationship between labile 
life-history traits. Testing the formal significance of the selec- 
tion differential can be performed by constraining the appro- 
priate covariance to zero and performing a likelihood ratio 
test between the unconstrained and constrained models. The 
procedure does not, however, allow description of nonlinear 
selection. In order to evaluate the potential for nonlinear se- 
lection on these properties, we suggest investigation of the 
map of relative fitness on the BLUP values for all reaction- 
norm properties. Under the paradigm of optimal reaction 
norms, one would expect that variation in slope maybe under 
stabilizing selection, which can hereby be graphically evalu- 
ated. 

All random regression models were solved using Restricted 
Maximum Likelihood (REML) in the programe ASReml 
(VSN International). This programe uses the delta method 
(see e.g., Lynch and Walsh 1998) to estimate the standard 
errors of functions of variances. The code used in this paper 
and example file of data is provided in the Supporting In- 
formation. In the implementation of the above equations in 
AsReml, the following need to be observed. (1) In random 
regression with a linear covariate, the variances of first-order 



and higher order terms are dependent on the scaling of the co- 
variate (e.g., Pinheiro and Bates 2000; Schaeffer 2004), which 
is dependent on the software used to implement the models. 
AsReml standardizes the covariate such that its minimal value 
is -1 and the average of the covariate values is zero. Hence, 
variances in elevation are defined for the average of all the 
covariate values (which may or may not equal the average 
covariate when taking the number of observations into ac- 
count). Random-regression slope is then defined as change 
in response variable per unit equal to the difference between 
minimum and mean value of the covariate. (2) The within- 
individual variance and variance in slope(s) of lifetime fitness 
cannot be constrained to exactly zero and must, in practice, 
be constrained to a very small value. See the Supporting In- 
formation for additional details. 

Example of the Random-Regression 
Approach: Ural Owl Clutch 
Size-Laying Date Reaction Norms 

Ural owl data 

A nest-box breeding population of Ural owls was studied in 
Paijat-Hame, southern Finland, in an area of 1500 km 2 dur- 
ing 1977-2009. Laying date of the first egg was established 
either by finding the incomplete clutch and calculating back 
assuming that eggs are laid at 2 days interval or by backdat- 
ing first the hatching date of the nestlings (based on a daily 
growth curve) and assuming an incubation period of 32 days 
for the first egg (largest nestling) . Any error in establishing 
laying date using these approaches is likely to be small in 
comparison to the length of the period of laying, which spans 
2 months. See Pietiainen (1989) and Kontiainen et al. (2010) 
for details. Each year, practically all Ural owls females were 
caught, either during incubation or when caring for their 
offspring. All females were aged by plumage characteristics 
to age categories of 1, 2, or >3 years old (Pietiainen and Kol- 
unen 1986) and were ringed with a metal ring in order to 
allow lifelong identification. All offspring were ringed and all 
boxes were checked after fledging to verify successful fledging 
of offspring. 

Data used in the analysis 

For the analyses, individuals that started their breeding career 
during the study period were included. We thus omitted in- 
dividuals trapped for the first time in 1977, but assumed that 
all females caught for the first time in one of our boxes from 
1978 onwards were first breeders. We included data on clutch 
size and laying date on all females that started their breeding 
career between 1978 and 2008 (31 years). We included age 
as a fixed-effect factor with three levels (1,2, 3+ years old). 
Laying date was standardized relative to the all-time median 
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laying date of 31 March (= 0), such that 1 April was 1 and 30 
March was -1. 

As an estimate of lifetime fitness, we used the total number 
of fledglings that a female has produced during her entire 
breeding career, which in this population correlates with the 
number of descendants after several generations (Brommer 
et al. 2004). Relative fitness w; of Equation (5) was calculated 
by dividing an individual's lifetime production of fledglings 
by the overall mean lifetime production of fledglings (fitness 
is relative to an average of unity) . 

We could analyze the clutch size-laying date relationships 
and their fitness consequences of 361 Ural owl females with 
1115 breeding attempts (3.1 records of clutch size and lay- 
ing date per female on average). All females with known 
laying date and clutch size that had bred at least once in our 
study population were included. Results did not change qual- 
itatively when animals with only one breeding record were 
excluded. 

Results: Interindividual Variation in 
the Clutch Size-Laying Date 
Relationship 

The annual mean Ural owl clutch size decreased approxi- 
mately linearly with advancing laying date (Fig. 2; Linear 
regression, intercept: 3.39 ± 0.098, f 30 = 34.5, P < 0.001; 
slope: -0.087 ± 0.014, t 30 = -6.3, P < 0.001; slope 2 : 0.00086 
± 0.0011, f 30 = 0.78, P = 0.44). Random regression analysis 
(Table 1 ) showed that there was substantial variation in clutch 



5 —i 




Laying date 

Figure 2. Yearly mean in clutch size and laying date for all Ural owl fe- 
males included in the analysis. Laying date was expressed in days relative 
to its long-term median (31 March). Data consist of 1 1 14 observations 
collected over 31 years. Line displays the linear regression (coefficients 
and test reported in the text). 



size across years that was not explained by variation in laying 
date or by age (model 2). Further, there was clear evidence of 
variance across females in all the elevation ind 0 terms (model 
3), linear slope terms ind\ (model 4), and quadratic slope 
terms indi (model 5) of their clutch size-laying date rela- 
tionships. The model that allowed individual-specific cubed 
[indi term) clutch size-laying date relationships (model 6, 
Table 1) was not a further improvement in model fit. We 
thus assumed that a second-order polynomial described the 
variation across individuals in their clutch size-laying date 
relationship most parsimoniously. Investigation of the fixed 
effects estimated under the final model (Table 1) showed that 
clutch size declined with 0.06 egg per day of advancing laying 
date and that 1 -year-old females tended to produce a smaller 
clutch size (although not significant, this effect was kept in 
the model). 

To visualize the thus-modeled relationship between clutch 
size and laying date, we plotted the variance in clutch size as 
a function of laying date (Fig. 3 A) based on the derivations 
derived by Fischer et al. (2004) . In addition, we plotted the in- 
dividual clutch size-laying date relationships for the females, 
based on the BLUP values from model 5 (Table 1) of ind 0 , 
indi, and ind 2 (Fig. 3B). Both methods showed clearly that 
there was considerable crossing of the reaction norms such 
that variance in clutch size across individuals was minimal 
around the average laying date (0) and high for both early 
and late clutches. 

Selection on the clutch size-laying date 
relationship 

Selection could be detected on all parameters that described 
the clutch size-laying date relationship (Table 2). Not sur- 
prisingly, the selection differential on the elevation of clutch 
size (expected clutch size at the average of laying date val- 
ues) showed clear selection for increased clutch size. There 
was an insignificant selection differential for a steeper (more 
negative) slope of the clutch size-laying date relationship. 
Lastly, there was significant positive selection for the second- 
order random regression term of the clutch size-laying date 
relationship. That is, female whose clutch size-laying date re- 
lationship was more concave had higher lifetime fitness. Be- 
cause these selection differentials take into account the whole 
(co)variance matrix (Table 2), selection on the second-order 
term occurred in addition to the selection on its correlate 
elevation. Because the bivariate random regression can only 
describe linear selection, we investigated the plots of rela- 
tive fitness versus the BLUP values for the reaction-norm 
properties to find evidence for stabilizing selection on linear 
slope or elevation. One may expect stabilizing selection on 
the reaction norm in case both shallower and steeper slopes 
reduce fitness (e.g., Fig. 1). Nevertheless, these plots revealed 
no evidence for stabilizing selection (Fig. SI). 
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Table 1. Hierarchical random regression linear mixed models of clutch size as a function of laying date of increasing polynomial order and their 
significance, as specified by Equation (1). Reported in the upper part of the table are the REML variances (covariances are reported in the text). The 
significance of higher order random regression is tested by a likelihood ratio test comparing the likelihood of each model to the hierarchical lower one. 
Residual variance is denoted for each model, and the coefficients of the random regression (Equation [1]) are listed in increasing order such that ind, 
is the first-order (linear) coefficient, ind 2 the second order. For each test, the degrees of freedom (df) are given by the number of additional variances 
and covariances that are estimated. The most parsimonious model is presented in bold. The lower part of the table presents the fixed effects and their 
standard error of the most parsimonious model, with their significance tested using f-tests. 



Random regression variances Test between models 



Model 


Residuals 


Year 


ind 0 


ind, 


ind 2 


ind 3 


LogL 


X 2 


df 


P 


1 


0.735 












-397.2 








2 


0.575 


0.179 










-295.2 


204.0 


1 


<0.001 


3 


0.452 


0.165 


0.126 








-264.8 


60.8 


1 


<0.001 


4 


0.413 


0.161 


0.122 


0.417 






-256.7 


16.2 


2 


0.003 


5 


0.397 


0.156 


0.141 


0.240 


0.975 




-251.7 


10.0 


3 


0.019 


6 


0.391 


0.156 


0.146 


0.234 


1.083 


0.487 


-249.3 


4.8 


4 


0.31 


Fixed effect 


Estimate 




df (nom) 


df(den) 


F 


P 










M 


3.38 ±0.080 




















Laying date 


-0.062 ±0.0032 




1 


382.3 


377.3 


<.001 










Age 


1: -0.22 ±0.14 
2: 0.032 ±0.10 




2 


983.6 


1.3 


0.285 











Interpretation of the findings 

We find that female Ural owls differ in the elevation, lin- 
ear and quadratic slopes of their clutch size-laying date re- 
lationships. This variation across individuals is important, 
because it creates the potential for selection to shape the 
plasticity in clutch size as a function of laying date. We in- 
deed show that selection operates on this plasticity, such that 
(1) females enjoy greater fitness when they have a larger el- 
evation. That is, a larger clutch size at the mean value of 
observed laying dates, which is a close correlate of an individ- 
ual's mean clutch size (cf. Brommer et al. 2003), (2) females 
whose clutch size declines steeper with advancing laying date 
and who show a stronger nonlinear (concave) curvature in 
their clutch size-laying date relationship have higher fitness. 
Inspection of the interindividual clutch size-laying date rela- 
tionship estimated by random regression (Fig. 3) shows that 
individuals with a concave curvature are those which have a 
relatively moderate clutch size at the early laying dates (before 
the median) and a relatively low clutch size during late laying 
dates (after the median). In the Ural owl, fledging success in 
brood sizes exceeding five offspring generally is poor, because 
hatching asynchrony increases rapidly with increasing clutch 
size (Kontiainen et al. 2010). For example, in a brood of 5, 
the smallest offspring is already 4-6 days younger than its 
oldest sibling and mortality of the youngest nestling in such 
a large brood is likely (Kontiainen et al. 2010). Therefore, the 
production of large clutch sizes adds little to lifetime fledgling 
production. In addition, late broods are typically found un- 
der poor environmental conditions (Pietiainen 1989) and it 
may be a good strategy to reproduce conservatively under 



those conditions. Thus, we believe that our selection anal- 
ysis captures biologically plausible processes in this species. 
The main finding is, nevertheless, that there is selection on 
properties that describe the plasticity of seasonal adjustment 
in clutch size. Provided part of this latter variation across fe- 
males in plasticity is heritable, selection has the potential to 
shape the clutch size-laying date relationship in this species. 

Discussion 

Theory states that the covariation in seasonal life-history 
traits such as clutch size and laying date arises from 
individual-specific fitness optima, suggesting that deviation 
from the reaction norm reduces fitness. Here, we advocate 
an approach based on random regression to first describe 
and test whether there is variation across individuals in the 
properties describing a trait-trait reaction norm and to fur- 
ther test for selection on these properties. We show that this 
approach indeed is effective in describing variation in Ural 
owl clutch size-laying date relationships and detects selection 
on the reaction-norm properties. Using a random regression 
mixed model to describe trait-trait reaction norms and se- 
lection on them avoids having to first extract information on 
reaction-norm properties from individual-specific regression 
and then regressing these estimates on fitness (e.g., Brommer 
et al. 2003; Nussey et al. 2005a). Hence, uncertainties in the 
estimate of an individual's reaction-norm properties are in- 
corporated and selection on the mean trait (i.e., elevation) is 
considered simultaneously to the within-individual covari- 
ation between traits as well as their codependency by con- 
sidering the complete matrix of variances and covariances 
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Figure 3. (A) Plot of the variance in clutch size (solid line) and its ap- 
proximate 95% confidence interval (dashed line) as a function of laying 
date. Values based on statistics of model 5 in Table 1 and applying the 
results of Fisher et al. (2004). (B) Plot of the reaction norms based on 
the fixed effects of elevation and slope and the Best Linear Unbiased 
Predictor (BLUP) values of model 5 for the individual females' reaction 
norms (Table 1). 



between the trait under investigation and fitness. The ap- 
proach works especially well on data consisting of repeated 
measures of traits on a set of individuals, of which the sea- 
sonal life-history traits in the case scenario considered is but 
one example. 

Selection acts on phenotypes, and, in case of phenotypic 
plasticity, the phenotype is an environment-specific realiza- 
tion of the genotype. Empirically measuring the fitness of 
the reaction norm requires integrating the fitness of all the 
phenotypes produced by a genotype over the entire environ- 
mental gradient (Weis and Gorman 1990). However, in the 
case of labile traits, we consider trait expression by a set of in- 
dividuals. Although each individual is a unique genotype, the 
variance between individuals in the properties describing the 
trait-trait reaction norms includes both the additive genetic 
variance and also other, nonheritable, consistent differences 
across individuals (so-called permanent environmental ef- 
fects, Lynch and Walsh 1998). In our example, environmental 
effect associated with individuals may arise because Ural owls 
are strongly territorial and essentially breed their entire career 
in the same territory (Saurola 1987). Thus, variation across 
territories is confounded with variation across individuals 
and it is therefore possible that the variation in reaction- 
norm properties we here document is entirely caused by 
environmental differences across individuals. Nevertheless, 
provided sufficient pedigree information is available, a ran- 
dom regression animal model can be used to further partition 
the phenotypic properties into their additive genetic versus 
permanent environmental components (Nussey et al. 2007). 
Such a random regression animal model would allow esti- 
mation of the additive genetic (co)variance matrix G for the 
reaction-norm properties. 

Nonlinear selection on reaction norm 
properties 

A major downside of the bivariate random regression ap- 
proach in quantifying selection on plasticity is that it does 
not allow for nonlinear selection. For a single trait, non- 
linear selection is estimated by calculating the covariance 
between the square of the trait and relative fitness (Falconer 
and MacKay 1996). The reaction-norm slope is an emergent 
property determined by the within -individual covariance be- 
tween the traits and cannot be directly measured. Selection 
on higher order polynomial terms, as we found in our ex- 
ample case for the second-order polynomial of the Ural owl 
clutch size-laying date relationship, refer to linear selection 
on specific curvature of the reaction norm itself rather than 
the shape of natural selection. Stabilizing selection on slope 
would be expected in case the reaction norm is maintained 
by selection. Nevertheless, there is generally a component 
of directional selection, except in the very restrictive case 
of completely symmetrical reductions in fitness if the slope 
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Table 2. Variance-covariance matrix (with SE) of the parameters involved in bivariate random regression of clutch size and lifetime fledgling production 
(LFP) as a function of laying date. LFP was expressed as relative fitness, by dividing it with its phenotypic mean (mean fitness of 1). The diagonal 
presents the variances in elevation (ind Q ), linear slope (ind,), and quadratic slope (ind 2 ), and LFP. Covariances are below the diagonal. The last row 
therefore presents the selection differential (S) on the reaction-norm properties. Correlations are reported above the diagonal. Model log-likelihood is 
-342.19. Significance of selection was tested using a likelihood ratio test between this full model and the model where the respective element of the 
matrix was constrained to zero. The log-likelihood (LogL) of the constrained model and the test statistics is reported in the bottom part of the table, 
together with the intensity if selection (i.e., 5 divided by the standard deviation of the respective parameter). The denominator degrees of freedom 
cannot be estimated numerically for this model and were conservatively equaled to the number of individuals. 





Bivariate random regression (co)variances 








ind 0 


ind, 


ind z 


LFP 


ind 0 


0.129 ±0.041 


-0.407 ± 0.23 


0.346 ±0.23 


0.394 ±0.11 




ind, 


-0.074 ± 0.043 


0.260 ±0.15 


0.218 ±0.33 


-0.194 ±0.16 






0.099 ±0.095 


0.089 ±0.13 


U.UtJ _!_ U.JJ 


\J . -) U _!_ \J . 1 O 




LFP 


0.104 ±0.031 


-0.073 ±0.056 


0.304 ± 0.092 


0.539 ±0.042 




Fixed effect on clutch 


size 


df (nom) 


df (den) 


F 


P 


M 


3.35 ± 0.079 










Laying date 


-0.060 ±0.0031 


1 


361 


379.4 


<0.001 


Age 


1: -0.19 ± 0.14 












2: 0.027 ±0.10 


2 


361 


0.88 


0.42 


Fixed effect on LFP 




df (nom) 


df (den) 


F 


P 


M 


1.54 ± 0.19 










Cohort year 




26 


361 


1.40 


0.095 


Parameter 


Intensity of selection (/) 


LRT 












LogL 


X 2 


df 


P 


indo 


0.29 


-340.6 


15.9 


1 


<0.001 


ind, 


-0.14 


-335.0 


1.83 


1 


0.17 


ind 2 


0.38 


-340.0 


13.4 


1 


<0.001 



deviates from the optimal reaction norm. Furthermore, com- 
parative studies underline that directional selection is typ- 
ically stronger than stabilizing selection (Kingsolver et al. 
2001; but see Stinchcombe et al. 2008). Investigation of the 
mapping of fitness to the BLUPs of the reaction-norm prop- 
erties, as we did here, should reveal whether stabilizing selec- 
tion on reaction norms may occur (although it does not allow 
a formal test). We advice that the estimates of selection on 
reaction-norm properties are considered with caution, but 
also emphasize that alternative approaches may lead to spu- 
rious findings, especially when the uncertainty of the slope 
estimate is ignored. 

Applications of the random-regression 
approach in evolutionary ecology 

The approach used here to quantify selection on trait-trait 
reaction norms is equally useful for reaction norms of trait 
versus environment. The typical use of random regression 
mixed models is to model the covariance of a trait and 
an abiotic variable. For example, reaction norms in laying 
date and temperature have been explored in several stud- 
ies in wild populations (e.g., Brommer et al. 2005, 2008; 
Nussey et al. 2005b; Charmantier et al. 2008; Husby et al. 



2010). Such random regression models can be readily ex- 
tended to a bivariate form where fitness is also consid- 
ered (Equation [6]) in order to analyze selection on the 
trait-environment reaction-norm properties. Understand- 
ing selection on trait-environment plasticity may be cru- 
cial, because under environmental change organisms may 
be selected to alter the plasticity in their trait-environment 
relationship, rather than merely their mean trait value. For 
example, Nussey et al. (2005b) showed that selection on the 
BLUPs of the slope of laying date-temperature relationships 
of great tits Varus major in the Netherlands increased as the 
local spring climate warmed during the last decades. 

Studies that consider repeated measures of different be- 
haviors collected on a set of individuals may benefit from the 
approach outlines in this paper. The level of integration of 
several behaviors within an individual is increasingly receiv- 
ing attention in ecology and evolutionary biology, because 
it captures "behavioural syndromes," consistency in behav- 
iors across different contexts (Sih et al. 2004). Under this 
paradigm, correlated behaviors (e.g., between "boldness" and 
"aggression") arise because individuals are relatively inflex- 
ible with respect to environmental conditions (Reale et al. 
2007). Hence, explicit consideration of "behavioural reac- 
tion norms" and understanding selection on their properties 
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has been recently advocated as a potentially fruitful line of 
research for understanding how natural selection shapes be- 
havioral syndromes (Dingemanse et al. 2010). In particular, 
the notion that behavioral syndromes are maintained in the 
population because plasticity, that is, reaction-norm slope, 
is genetically constrained (Sih et al. 2004) as opposed to be- 
ing shaped by selection can be explicitly explored using the 
random regression approach advocated here. 

Conclusion 

Phenotypic plasticity is a central concept in evolutionary ecol- 
ogy (Schlichting 1986; Scheiner 1993; Schlichting and Pigli- 
ucci 1998). The seminal work of Weis and Gorman (1990) 
has clarified that quantification of natural selection on the 
level of the reaction norm is pivotal in understanding the 
evolution of plasticity. From this perspective, repeated ob- 
servations on long-lived organisms under natural conditions 
in combination with the methodology we here presented 
constitute a powerful approach to quantify variation in the 
reaction-norm properties and linking them to variation in 
fitness. 
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